#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include "toroidal.h"

double innerproduct(fiq a, fiq b)
{
    double* ia = a.get_fi();
    double* ib = b.get_fi();
    double* qa = a.get_fq();
    double* qb = b.get_fq();
    
    double ret = 0.0;
    for(int i = 0; i < Nx*Ny; i++)
    {
        double rloc = (i/Ny)*dr + Rmin;
        ret += rloc*(ia[i]*ib[i] + qa[i]*qb[i])*dr*dz;
    }
    return ret;
}

double inner_vf(vfiq F1, vfiq F2)
{
    double* ifr1 = F1.get_ifr();
    double* ifz1 = F1.get_ifz();
    double* qfr1 = F1.get_qfr();
    double* qfz1 = F1.get_qfz();
    double* ifr2 = F2.get_ifr();
    double* ifz2 = F2.get_ifz();
    double* qfr2 = F2.get_qfr();
    double* qfz2 = F2.get_qfz();
    
    double ret = 0.0;
    for(int i = 0; i < Nx*Ny; i++)
    {
        double rloc = (i/Ny)*dr + Rmin;
        ret += rloc*(ifr1[i]*ifr2[i] + qfr1[i]*qfr2[i])*dr*dz;
        ret += rloc*(ifz1[i]*ifz2[i] + qfz1[i]*qfz2[i])*dr*dz;
    }
    return ret;
}

void eigen(double* data, int dim)
{
    gsl_matrix_view m = gsl_matrix_view_array(data, dim, dim);
    gsl_vector *eval = gsl_vector_alloc(dim);
    gsl_matrix *evec = gsl_matrix_alloc(dim, dim);
    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(dim);

    gsl_eigen_symmv(&m.matrix, eval, evec, w);
    gsl_eigen_symmv_free(w);
    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);

    for(int i = 0; i < dim; i++)
    {
        double eval_i = gsl_vector_get(eval, i);
        gsl_vector_view evec_i = gsl_matrix_column(evec, i);
        printf ("eigenvalue = %g\n", eval_i);
        printf ("eigenvector = \n");
    }
    gsl_vector_free(eval);
    gsl_matrix_free(evec);
}

void eigen_gen(double* A, double* B, int dim, int eign, double* coefs)
{
    gsl_matrix_view mA = gsl_matrix_view_array(A, dim, dim);
    gsl_matrix_view mB = gsl_matrix_view_array(B, dim, dim);
    gsl_vector *eval = gsl_vector_alloc(dim);
    gsl_matrix *evec = gsl_matrix_alloc(dim, dim);
    gsl_eigen_gensymmv_workspace *w = gsl_eigen_gensymmv_alloc(dim);

    gsl_eigen_gensymmv(&mA.matrix, &mB.matrix, eval, evec, w);
    gsl_eigen_gensymmv_free(w);
    gsl_eigen_gensymmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);

    for(int i = 0; i < dim; i++)
    {
        double eval_i = gsl_vector_get(eval, i);
        gsl_vector_view evec_i = gsl_matrix_column(evec, i);
        if(i>4)
            continue;
        printf ("\neigenvalue = %g\n", eval_i);
        for(int j = 0; j < dim; j++)
            printf("%f, ", gsl_vector_get(&evec_i.vector, j));
    }
    printf("\n");
    gsl_vector_view evec_get = gsl_matrix_column(evec, eign);
    for(int i = 0; i < dim; i++)
        coefs[i] = gsl_vector_get(&evec_get.vector, i);
    gsl_vector_free(eval);
    gsl_matrix_free(evec);
}

void linear_suppos(double* coefs, fiq* eigvecs, int n, fiq* eigm)
{
    double* outi = eigm->get_fi();
    double* outq = eigm->get_fq();
    for(int i = 0; i < n; i++)
    {
        double* ini = eigvecs[i].get_fi();
        double* inq = eigvecs[i].get_fq();
        for(int j = 0; j < Nx*Ny; j++)
        {
            outi[j] += coefs[i]*ini[j];
            outq[j] += coefs[i]*inq[j];
        }
    }
}